Ephedra regional gradient (ERG) analyses

ecoblender

alex.filazzola

Abstract

An under-examined component of the shrub-annual relationship is how regional drivers, such as climate, may alter the sign or magnitude of positive interactions. Interspecific interactions between plants have been shown to be strongly linked to climate, particularly temperature and precipitation. The stress-gradient hypothesis (SGH) predicts that higher abiotic stress (i.e. for deserts, increasing temperature and reduced precipitation) will increase the frequency of positive interactions among shrubs and their annual understory. Regional climate gradients also have indirect effects on plant composition, such as determining consumer abundance and soil nutrient composition. Nutrient availability is particularly affected by precipitation because of altered decomposition rates of organic matter and mineralization. Therefore, the strength of facilitation and operating mechanism of a shrub on the annual plant community may change along a regional gradient.

Hypothesis

We tested the hypothesis that positive interactions among shrubs and annual plants will increase with abiotic stress and reduce nutrient availability along a regional gradient of aridity.

Methods

Weather for growing season

Climate patterns within study

season1.sjd <- season1 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season1.mnp <- season1 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.sjd <- season2 %>% filter(Gradient<4) %>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
season2.mnp <- season2 %>% filter(Gradient>3)%>% group_by(year, month,days) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))

## Rain vs Temperature in 2016
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot1 <- barplot(height=season1.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot1[,1], season1.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot1 <- barplot(height=season1.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
points(plot1[,1], season1.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature at all sites (°C)", 4, line=3)

## Rain vs Temperature in 2017
par(mfrow=c(2,1))
par(mar=c(3.5,4.5,1,4.5))
plot2 <- barplot(height=season2.sjd$Precip, ylim=c(0,14), ylab="Average precipitation San Joaquin (cm)")
points(plot2[,1], season2.sjd$min.temp, type="l", col="#FF000050", lwd=2)
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature San Joaquin (°C)", 4, line=3)
par(mar=c(4.5,4.5,0,4.5))
plot2 <- barplot(height=season2.mnp$Precip, ylim=c(0,14), ylab="Average precipitation Mojave (cm)")
points(plot2[,1], season2.mnp$min.temp, type="l", col="#FF000050", lwd=2)
axis(1, plot1[c(1,30,60,90,120,150,180)], c("Nov","Dec","Jan","Feb","Mar","Apr","May"))
axis(4, at=seq(0,14,2), lab=seq(0,14,2), ylab="")
mtext("Average temperature Mojave (°C)", 4, line=3)

### 2016 The rain was inconsistent and mostly absent in the Mojave. This resulted in low germination and producitivty at the southern sites
### 2017 The rain was more plentiful, but in the northern sites, there appears to be a frost period after the majority of the rainfall. Need to check number of frost days

season1.frost <- season1 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season1.frost)
##              Site frost.days
## 1         Barstow  24.725275
## 2          Cuyama  29.670330
## 3   HeartofMojave  25.824176
## 4    PanocheHills  10.439560
## 5 SheepholeValley   6.043956
## 6          Tecopa  23.626374
## 7      TejonRanch  50.549451
season2.frost <- season2 %>% group_by(Site) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100)
data.frame(season2.frost)
##              Site frost.days
## 1         Barstow  17.679558
## 2          Cuyama  19.889503
## 3   HeartofMojave  16.574586
## 4    PanocheHills  14.917127
## 5 SheepholeValley   1.785714
## 6          Tecopa  25.966851
## 7      TejonRanch  35.911602
## Both years had comparable number of frost days

## Compare number of consecutive frost days (i.e. frost periods)
season1[,"frost"] <- ifelse(season1$min.temp<0, -99,season1$min.temp) ## identified days below freezing
season2[,"frost"] <- ifelse(season2$min.temp<0, -99,season2$min.temp) ## identified days below freezing
count.consec <- function(x) {max(rle(as.character(x))$lengths)}

season1.frost <- season1 %>% group_by(Site)  %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                  10
## 2          Cuyama                  10
## 3   HeartofMojave                  12
## 4    PanocheHills                   5
## 5 SheepholeValley                   7
## 6          Tecopa                  11
## 7      TejonRanch                  14
season2.frost <- season2 %>% group_by(Site) %>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   7
## 4    PanocheHills                   6
## 5 SheepholeValley                   3
## 6          Tecopa                   7
## 7      TejonRanch                  13
## compare only after plants have germinated
season1.frost <- season1 %>% group_by(Site) %>% filter(year>2015) %>% summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100, avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season1.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  12.396694     5.750413
## 2          Cuyama  11.570248     3.352893
## 3   HeartofMojave  16.528926     4.542149
## 4    PanocheHills   2.479339     6.777686
## 5 SheepholeValley   2.479339     9.394167
## 6          Tecopa  11.570248     6.342149
## 7      TejonRanch  33.884298     1.438017
season2.frost <- season2 %>% group_by(Site) %>%  filter(year>2016) %>%  summarize(frost.days=sum(min.temp<0, na.rm=T)/length(min.temp)*100,avg.min.temp=mean(min.temp, na.rm=T))
data.frame(season2.frost)
##              Site frost.days avg.min.temp
## 1         Barstow  11.666667     6.052500
## 2          Cuyama  16.666667     3.624167
## 3   HeartofMojave   8.333333     5.637838
## 4    PanocheHills   8.333333     6.065833
## 5 SheepholeValley   0.000000     9.386916
## 6          Tecopa  20.833333     5.938333
## 7      TejonRanch  28.333333     2.535833
season1.frost <- season1 %>% group_by(Site)  %>% filter(year>2015) %>% summarize(count.consec(frost))
data.frame(season1.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   4
## 4    PanocheHills                   2
## 5 SheepholeValley                   2
## 6          Tecopa                   4
## 7      TejonRanch                  10
season2.frost <- season2 %>% group_by(Site)  %>%  filter(year>2016)%>% summarize(count.consec(frost))
data.frame(season2.frost)
##              Site count.consec.frost.
## 1         Barstow                   5
## 2          Cuyama                   6
## 3   HeartofMojave                   3
## 4    PanocheHills                   3
## 5 SheepholeValley                   2
## 6          Tecopa                   5
## 7      TejonRanch                  10

Microenvironmental differences

We compared temperature and relative humidity between shrub and open microsites among all sites along the regional gradient


    Shapiro-Wilk normality test

data:  fit1$residuals
W = 0.87881, p-value < 2.2e-16
                     Df Sum Sq Mean Sq F value   Pr(>F)    
Microsite             1   26.5   26.51  59.886 1.24e-14 ***
gradient              6  233.9   38.98  88.065  < 2e-16 ***
Microsite:gradient    6   13.1    2.18   4.927 4.95e-05 ***
Residuals          4290 1898.9    0.44                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30 observations deleted due to missingness
Analysis of Deviance Table

Model: binomial, link: logit

Response: RH/100

Terms added sequentially (first to last)

                   Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
NULL                                4333     974.99             
Microsite           1     1.33      4332     973.66   0.2492    
gradient            6   528.18      4326     445.48   <2e-16 ***
Microsite:gradient  6     7.69      4320     437.79   0.2613    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Deviance Table

Model: Gamma, link: inverse

Response: swc

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev        F Pr(>F)    
NULL                             404    231.405                    
Site            6  189.683       398     41.722 381.2179 <2e-16 ***
Microsite       1    0.287       397     41.435   3.4579 0.0637 .  
Site:Microsite  6    0.696       391     40.739   1.3994 0.2136    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Deviance Table

Model: Gamma, link: inverse

Response: swc

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev        F Pr(>F)    
NULL                             419    183.193                    
Site            6  144.930       413     38.263 279.4537 <2e-16 ***
Microsite       1    0.080       412     38.183   0.9284 0.3358    
Site:Microsite  6    0.535       406     37.647   1.0322 0.4037    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

site level means

## site coordintes
site.gps <- read.csv("Data/ERGsites.csv")
# nutrient data for each site
nutrients <- read.csv("Data/ERG.soilnutrients.csv")
 ## long term climate data extracted from worldclim
clim.dat <- read.csv("Data/ERG.worldclim.csv")


## check correlation among long-term climate variables
cor(clim.dat[,3:8]) ## Temp.wettest.QR least correlated 
##                    Annual.Temp Temp.Seasonality Temp.wettest.QR
## Annual.Temp          1.0000000        0.9051791       0.7260617
## Temp.Seasonality     0.9051791        1.0000000       0.4650548
## Temp.wettest.QR      0.7260617        0.4650548       1.0000000
## Annual.precip       -0.8911069       -0.9650533      -0.5313460
## Precip.seasonality  -0.9715419       -0.9774646      -0.6165170
## Precipt.wettest.QR  -0.8839124       -0.9570957      -0.5523347
##                    Annual.precip Precip.seasonality Precipt.wettest.QR
## Annual.Temp           -0.8911069         -0.9715419         -0.8839124
## Temp.Seasonality      -0.9650533         -0.9774646         -0.9570957
## Temp.wettest.QR       -0.5313460         -0.6165170         -0.5523347
## Annual.precip          1.0000000          0.9613729          0.9980959
## Precip.seasonality     0.9613729          1.0000000          0.9543657
## Precipt.wettest.QR     0.9980959          0.9543657          1.0000000
## Obtain mean shrub traits for each site
shrubs <- read.csv("Data/ERG.shrub.csv")
shrubs <- subset(shrubs, Microsite=="shrub")
shrubs.mean <- shrubs %>% group_by(Gradient,Site) %>% summarise_all(funs(mean))
shrubs.mean <- data.frame(shrubs.mean)
shrubs.vars <- shrubs.mean[,c("Site","Gradient","volume","canopy","Dx","DxEph","Compaction")] 

## test correlation among shrub traits
cor(shrubs.vars[,3:7]) ## compaction x volume & DxEph x Dx high correlation
##                volume     canopy          Dx      DxEph  Compaction
## volume      1.0000000  0.2243902  0.40036651  0.2631563 -0.64927168
## canopy      0.2243902  1.0000000 -0.68479570 -0.3538094 -0.20833287
## Dx          0.4003665 -0.6847957  1.00000000  0.7803527 -0.06531485
## DxEph       0.2631563 -0.3538094  0.78035272  1.0000000  0.33769574
## Compaction -0.6492717 -0.2083329 -0.06531485  0.3376957  1.00000000
vifstep(shrubs.vars[,3:7], th=10) ## Dx showing collinearity problems
## 1 variables from the 5 input variables have collinearity problem: 
##  
## Dx 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( Compaction ~ canopy ):  -0.2083329 
## max correlation ( Compaction ~ volume ):  -0.6492717 
## 
## ---------- VIFs of the remained variables -------- 
##    Variables      VIF
## 1     volume 4.043473
## 2     canopy 1.474526
## 3      DxEph 2.874621
## 4 Compaction 3.760284
shrub.site <- shrubs.vars[,-c(1,2,5)]
rownames(shrub.site) <- shrubs.vars[,1]

## PCA of shrub characteristics
pca1 <- prcomp(log(shrub.site), scale=T)
plot(pca1)

biplot(pca1, scale = T)

summary(pca1) ## 77% variation explained
## Importance of components:
##                           PC1    PC2    PC3     PC4
## Standard deviation     1.2958 1.1879 0.8922 0.33722
## Proportion of Variance 0.4198 0.3528 0.1990 0.02843
## Cumulative Proportion  0.4198 0.7725 0.9716 1.00000
## PCA of site characteristics for season1
## Obtain mean weather variables for each site
season1.mean <- season1 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip),wind=mean(wind.speed, na.rm=T))
season1.mean  <- data.frame(season1.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress

##combine nutrients and long-term averages
site.vars <- data.frame(season1.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations, 
row.names(site.vars) <- shrubs.vars[,1]
cor(site.vars)
##              temp.var     Precip        wind  elevation
## temp.var   1.00000000 -0.8679054  0.05816084 -0.2699588
## Precip    -0.86790537  1.0000000 -0.37279135  0.1222504
## wind       0.05816084 -0.3727914  1.00000000  0.1054497
## elevation -0.26995879  0.1222504  0.10544974  1.0000000
site.vars2016 <- site.vars

##  check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum
## No variable from the 4 input variables has collinearity problem. 
## 
## The linear correlation coefficients ranges between: 
## min correlation ( wind ~ temp.var ):  0.05816084 
## max correlation ( Precip ~ temp.var ):  -0.8679054 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1  temp.var 6.640710
## 2    Precip 7.316232
## 3      wind 1.739558
## 4 elevation 1.142681
pca1 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca1)

biplot(pca1)

## check contribution of loadings
pca1$rotation
##                  PC1        PC2         PC3        PC4
## temp.var   0.6001002 -0.1669018  0.50554765 -0.5970303
## Precip    -0.6703313 -0.1105034 -0.09936844 -0.7270288
## wind       0.2743589  0.7981791 -0.43410209 -0.3149488
## elevation -0.3395039  0.5681927  0.73898773  0.1256633
aload <- abs(pca1$rotation)
sweep(aload, 2, colSums(aload), "/")
##                 PC1       PC2        PC3        PC4
## temp.var  0.3184748 0.1015355 0.28433406 0.33832382
## Precip    0.3557466 0.0672253 0.05588758 0.41199111
## wind      0.1456030 0.4855763 0.24415109 0.17847449
## elevation 0.1801756 0.3456629 0.41562726 0.07121058
summary(pca1) ## 89% variation explained
## Importance of components:
##                          PC1    PC2    PC3     PC4
## Standard deviation     1.456 1.0458 0.8577 0.22479
## Proportion of Variance 0.530 0.2734 0.1839 0.01263
## Cumulative Proportion  0.530 0.8034 0.9874 1.00000
gradient1.season1 <- pca1$x[,1]
gradient2.season1 <- pca1$x[,2]

## season2
season2.mean <- season2 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip, na.rm=T),wind=mean(wind.speed, na.rm=T))
season2.mean  <- data.frame(season2.mean)
## extract key variables
## dropped RH, and chose min.temp because least correlation with others & cold stress

##combine nutrients and long-term averages
site.vars <- data.frame(season2.mean[,3:5],site.gps["elevation"]) ## drop Phosphorus and other climate variables because of correlations, 
row.names(site.vars) <-  shrubs.vars[,1]
cor(site.vars)
##             temp.var     Precip       wind  elevation
## temp.var   1.0000000 -0.7427028 -0.2262409 -0.3715471
## Precip    -0.7427028  1.0000000 -0.1551517  0.1813600
## wind      -0.2262409 -0.1551517  1.0000000  0.1052805
## elevation -0.3715471  0.1813600  0.1052805  1.0000000
site.vars2017 <- site.vars


##  check for collinearity
vifstep(site.vars, th=10) ## remove potassium and temperature minimum
## No variable from the 4 input variables has collinearity problem. 
## 
## The linear correlation coefficients ranges between: 
## min correlation ( elevation ~ wind ):  0.1052805 
## max correlation ( Precip ~ temp.var ):  -0.7427028 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables      VIF
## 1  temp.var 3.438786
## 2    Precip 3.035337
## 3      wind 1.402010
## 4 elevation 1.192010
pca2 <- prcomp(log(abs(site.vars)), scale=T)
plot(pca2)

biplot(pca2)

summary(pca2) ## 85% variation explained
## Importance of components:
##                          PC1    PC2    PC3     PC4
## Standard deviation     1.434 1.0121 0.8809 0.37883
## Proportion of Variance 0.514 0.2561 0.1940 0.03588
## Cumulative Proportion  0.514 0.7701 0.9641 1.00000
## check contribution of loadings
pca2$rotation
##                  PC1         PC2         PC3       PC4
## temp.var   0.6557532 -0.09123357  0.20565681 0.7206729
## Precip    -0.6176537 -0.16891771 -0.40105332 0.6550778
## wind      -0.1019884  0.97279862  0.06826009 0.1964733
## elevation -0.4220070 -0.12963831  0.89005734 0.1135866
aload <- abs(pca2$rotation)
sweep(aload, 2, colSums(aload), "/")
##                  PC1        PC2       PC3        PC4
## temp.var  0.36483385 0.06695609 0.1314078 0.42749339
## Precip    0.34363685 0.12396828 0.2562596 0.38858328
## wind      0.05674213 0.71393442 0.0436159 0.11654531
## elevation 0.23478717 0.09514122 0.5687167 0.06737801
## define gradients
gradient1.season2 <- pca2$x[,1]
gradient2.season2 <- pca2$x[,2]

Compare gradient to phytometer

# crs.world <-CRS("+proj=longlat +datum=WGS84")
# gps <- data.frame(x=site.gps$long,y=site.gps$lat)
# coordinates(gps) <- ~x+y
# proj4string(gps) <- crs.world
# 
# ## Download and extract PET/aridity
# r1 <- raster("C:\\Users\\Alessandro\\Downloads\\Global Aridity - Annual\\AI_annual\\ai_yr\\w001001.adf", package="raster") #aridity
# r2 <- raster("C:\\Users\\Alessandro\\Downloads\\Global PET - Annual\\PET_he_annual\\pet_he_yr\\w001001.adf", package="raster") #potential evapotranspiration
# calclim <- stack(r1,r2)
# names(calclim) <- c("aridity","PET")
# arid.vals <- raster::extract(calclim, gps)
# rownames(arid.vals) <- site.gps[,"name"]
# colnames(arid.vals)[1] <- "Gradient"
# write.csv(arid.vals, "Data/aridity.PET.csv")
arid.vals <- read.csv("Data/aridity.PET.csv")
colnames(arid.vals)[1] <- "Site"


par(mfrow=c(1,2))
mean.phyto <- census %>% filter(Census=="end") %>%group_by(Year, Gradient, Site, Microsite) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
mean.phyto <- data.frame(mean.phyto)

rii.data <- rii(census, 1:6, var=9:15)
rii.mean <- rii.data %>% group_by(Year, Gradient,Site) %>% summarise_if(is.numeric, funs(mean(., na.rm=T)))
rii.mean <- data.frame(rii.mean)

season1.mean <- season1 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip),wind=mean(wind.speed, na.rm=T), max.temp=abs(mean(max.temp, na.rm=T)),min.temp=abs(mean(min.temp, na.rm=T)),avg.temp=abs(mean(avg.temp, na.rm=T)))
s1.mean <- data.frame(season1.mean)

season2.mean <- season2 %>% group_by(Gradient,Site) %>% summarise(temp.var=var(avg.temp,na.rm=T),Precip=sum(Precip, na.rm=T),wind=mean(wind.speed, na.rm=T), max.temp=abs(mean(max.temp, na.rm=T)),min.temp=abs(mean(min.temp, na.rm=T)),avg.temp=abs(mean(avg.temp, na.rm=T)))
s2.mean <- data.frame(season2.mean)


s1.mean[,"arid.gradient"] <- log(s1.mean[,"Precip"]/arid.vals[,"PET"])
s2.mean[,"arid.gradient"] <- log(s2.mean[,"Precip"]/(arid.vals[,"PET"]))

## collect aridity gradient data
site.climate <- rbind(s1.mean,s2.mean)
site.climate[,"Year"] <- as.factor(c(rep("2016",7),rep("2017",7)))

## combine climate data with census
mean.phyto  <- merge(mean.phyto , site.climate, by=c("Year","Site"))

phyto.2016 <- subset(mean.phyto, Year==2016)

plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"],1), phyto.2016[phyto.2016$Microsite=="open","phyto.biomass"], ylim=c(0,2), ylab=c("all phytometer biomass (g)"), xlab=c("aridity gradient"))
points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"],2), phyto.2016[phyto.2016$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,1.8, c("shrub","open"), pch=c(19,1))
text(-2.3, 1.8, "2016", cex=1.5)

plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"], rii.mean[rii.mean$Year==2016,"phyto.biomass"], ylim=c(-0.4,0.4), ylab="RII all phytometer biomass", pch=19, xlab=c("aridity gradient"))
text(-2.3, 0.38, "2016", cex=1.5)

phyto.2017 <- subset(mean.phyto, Year==2017)

plot(jitter(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"],1), phyto.2017[phyto.2017$Microsite=="open","phyto.biomass"], ylim=c(0,3), xlab=c("aridity gradient"),ylab=c("all phytometer biomass (g)"))
points(jitter(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"],2), phyto.2017[phyto.2017$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,1.8, c("shrub","open"), pch=c(19,1))
text(-2.0, 2.8, "2017", cex=1.5)

plot(phyto.2017[phyto.2017$Microsite=="open","arid.gradient"], rii.mean[rii.mean$Year==2017,"phyto.biomass"], ylim=c(-0.4,0.4), ylab="RII all phytometer biomass", pch=19, xlab=c("aridity gradient"))
text(-2.0, 0.38, "2017", cex=1.5)

## regress gradient and microsite on biomass for phytometers in 2016
m1 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m1) ## nothing significant
## 
## Call:
## lm(formula = ihs(Phacelia.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26646 -0.16692  0.01821  0.09053  0.43991 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                   0.75896    0.31899   2.379   0.0387 *
## arid.gradient                 0.16726    0.09147   1.829   0.0974 .
## Micrositeshrub                0.81930    0.45112   1.816   0.0994 .
## arid.gradient:Micrositeshrub  0.19010    0.12936   1.470   0.1724  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2383 on 10 degrees of freedom
## Multiple R-squared:  0.674,  Adjusted R-squared:  0.5762 
## F-statistic: 6.893 on 3 and 10 DF,  p-value: 0.008498
m2 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m2) ## nothing significant
## 
## Call:
## lm(formula = ihs(Plantago.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2016)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.070995 -0.048467  0.003316  0.024049  0.115472 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                   0.27085    0.07901   3.428  0.00646 **
## arid.gradient                 0.05886    0.02266   2.598  0.02658 * 
## Micrositeshrub               -0.12531    0.11174  -1.121  0.28831   
## arid.gradient:Micrositeshrub -0.02928    0.03204  -0.914  0.38239   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05904 on 10 degrees of freedom
## Multiple R-squared:  0.4794, Adjusted R-squared:  0.3232 
## F-statistic: 3.069 on 3 and 10 DF,  p-value: 0.07774
m3 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2016)
summary(m3)## aridity significant
## 
## Call:
## lm(formula = ihs(Salvia.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15051 -0.09099 -0.01568  0.03332  0.26045 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                   0.748316   0.194293   3.851   0.0032 **
## arid.gradient                 0.159284   0.055714   2.859   0.0170 * 
## Micrositeshrub               -0.000770   0.274771  -0.003   0.9978   
## arid.gradient:Micrositeshrub  0.004763   0.078792   0.060   0.9530   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1452 on 10 degrees of freedom
## Multiple R-squared:  0.6281, Adjusted R-squared:  0.5165 
## F-statistic:  5.63 on 3 and 10 DF,  p-value: 0.01599
## regress gradient and microsite on biomass for phytometers in 2017
m4 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m4) ## nothing significant
## 
## Call:
## lm(formula = ihs(Phacelia.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2017)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53579 -0.12720  0.00415  0.29889  0.40793 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                    -6.725      3.895  -1.727    0.135
## arid.gradient                  -2.514      1.315  -1.912    0.104
## Micrositeshrub                  6.929      4.012   1.727    0.135
## arid.gradient:Micrositeshrub    2.341      1.363   1.718    0.137
## 
## Residual standard error: 0.4168 on 6 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.3953, Adjusted R-squared:  0.09289 
## F-statistic: 1.307 on 3 and 6 DF,  p-value: 0.3557
m5 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m5) ## aridity significant
## 
## Call:
## lm(formula = ihs(Plantago.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2017)
## 
## Residuals:
##        15        16        19        20        23        24        25 
## -0.008390 -0.070144  0.006678 -0.031325  0.126617  0.249998 -0.124905 
##        26 
## -0.148529 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   -1.4847     1.6331  -0.909    0.415
## arid.gradient                 -0.6486     0.5513  -1.177    0.305
## Micrositeshrub                 3.1150     2.3096   1.349    0.249
## arid.gradient:Micrositeshrub   1.1087     0.7796   1.422    0.228
## 
## Residual standard error: 0.1748 on 4 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.4909, Adjusted R-squared:  0.109 
## F-statistic: 1.285 on 3 and 4 DF,  p-value: 0.3935
m6 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2017)
summary(m6)## nothing significant
## 
## Call:
## lm(formula = ihs(Salvia.biomass) ~ arid.gradient * Microsite, 
##     data = phyto.2017)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52276 -0.13176 -0.02786  0.05787  0.53706 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)
## (Intercept)                   -0.5442     0.6266  -0.869    0.405
## arid.gradient                 -0.4200     0.2423  -1.733    0.114
## Micrositeshrub                 0.2852     0.8861   0.322    0.754
## arid.gradient:Micrositeshrub   0.1488     0.3427   0.434    0.673
## 
## Residual standard error: 0.335 on 10 degrees of freedom
## Multiple R-squared:  0.3112, Adjusted R-squared:  0.1046 
## F-statistic: 1.506 on 3 and 10 DF,  p-value: 0.2723
# 
# ## Phacelia season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="shrub","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Phacelia.biomass"], pch=19)
# 
# m1 <- lm(ihs(Phacelia.biomass)~arid.gradient * Microsite, data=phyto.2016)
# m2 <- lm(ihs(Plantago.biomass)~arid.gradient * Microsite, data=phyto.2016)
# m3 <- lm(ihs(Salvia.biomass)~arid.gradient * Microsite, data=phyto.2016)
# 
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Phacelia.biomass"], ylim=c(-0.4,0.4))
# 
# ## Plantago season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Plantago.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Plantago.biomass"], pch=19)
# 
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Plantago.biomass"], ylim=c(-0.4,0.4))
# 
# ## Salvia season 1
# plot(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,1), phyto.2016[phyto.2016$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] ,2), phyto.2016[phyto.2016$Microsite=="shrub","Salvia.biomass"], pch=19)
# 
# plot(phyto.2016[phyto.2016$Microsite=="open","arid.gradient"] , rii.mean[rii.mean$Year==2016,"Salvia.biomass"], ylim=c(-0.4,0.4))
# 
# 
# 
# ## Phacelia season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Phacelia.biomass"], pch=19)
# 
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Phacelia.biomass"], ylim=c(-0.4,0.4))
# 
# 
# ## Plantago season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Plantago.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Plantago.biomass"], pch=19)
# 
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Plantago.biomass"], ylim=c(-0.4,0.4))
# 
# ## Salvia season 2
# plot(jitter(gradient1.season2,1), phyto.2017[phyto.2017$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(gradient1.season2,2), phyto.2017[phyto.2017$Microsite=="shrub","Salvia.biomass"], pch=19)
# 
# 
# plot(gradient1.season2, rii.mean[rii.mean$Year==2017,"Salvia.biomass"], ylim=c(-0.4,0.4))


## both seasons


## responses
plot(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"], mean.phyto[mean.phyto$Microsite=="open","phyto.biomass"], ylim=c(0,3), ylab=c("all phytometer biomass (g)"))
points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="shrub","phyto.biomass"], pch=19)
legend(-4.5,2.8, c("shrub","open"), pch=c(19,1))
text(-2.3, 2.8, "both years", cex=1.3)

# 
# ## Phacelia
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Phacelia.biomass"], ylim=c(0,2))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Phacelia.biomass"], pch=19)
# m1 <- lm(ihs(Phacelia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m1)
# 
# ## Plantago 
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Plantago.biomass"], ylim=c(0,0.5))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Plantago.biomass"], pch=19)
# m2 <- lm(ihs(Phacelia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m2)
# 
# 
# ## Salvia season 1
# plot(jitter(mean.phyto[mean.phyto$Microsite=="open","arid.gradient"],1), mean.phyto[mean.phyto$Microsite=="open","Salvia.biomass"], ylim=c(0,2))
# points(jitter(mean.phyto[mean.phyto$Microsite=="shrub","arid.gradient"],2), mean.phyto[mean.phyto$Microsite=="shrub","Salvia.biomass"], pch=19)
# m3 <- lm(ihs(Salvia.biomass) ~ arid.gradient * Microsite, data=mean.phyto)
# summary(m3)



### Rii plots

rii.mean <- merge(rii.mean, site.climate, by=c("Year","Site"))


## biomass by RII
plot(mean.phyto[mean.phyto$Microsite=="open","Phacelia.biomass"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.3,0.3), pch=19, ylab="RII biomass", xlab="average phytometer biomass per site")
points(mean.phyto[mean.phyto$Microsite=="open","Plantago.biomass"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
points(mean.phyto[mean.phyto$Microsite=="open","Salvia.biomass"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
abline(h=0, lty=2, lwd=2)
legend(0.38,-0.15, c("Phacelia","Plantago","Salvia"), pch=21, pt.bg=c("black","blue","red"))

## rii variability


plot(rii.mean[,"arid.gradient"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.5,0.5), pch=19, ylab="Rii phytometer biomass", xlab="aridity gradient")
points(rii.mean[,"arid.gradient"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
points(rii.mean[,"arid.gradient"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
abline(h=0, lty=2, lwd=2)
legend(-4.5, 0.45, c("Phacelia","Plantago","Salvia"), pch=21, pt.bg=c("black","blue","red"))

## compare convex hull
## hull function
Plot_ConvexHull<-function(xcoord, ycoord, lcolor){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor)
}  

Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Phacelia.biomass"], lcolor = "black")
Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Plantago.biomass"], lcolor = "blue")
Plot_ConvexHull(xcoord = rii.mean[,"arid.gradient"], ycoord = rii.mean[,"Salvia.biomass"], lcolor = "red")


r.df <- rii.mean[,c("arid.gradient","Phacelia.biomass","Plantago.biomass","Salvia.biomass")]


## transform data into long format
r.df <- gather(r.df, species, biomass, 2:4)


## function to create a convex hull around data
hull.time <- function(df){
  df[chull(df$arid.gradient,df$biomass),]  # chull really is useful, even outside of contrived examples.
  }

## break data into lists and recombined to apply polygon
splitData <- split(r.df, r.df$species)
appliedData <- lapply(splitData, hull.time)
combinedData <- do.call(rbind, appliedData)

## create polygon 
poly1 <- ggplot(r.df, aes(x=arid.gradient, y=biomass))+
geom_polygon(data = combinedData,  # This is also a nice example of how to plot
                          aes(x=arid.gradient, y=biomass, group = species, color=species), fill=NA,  # two superimposed geoms
                          alpha = 1/2)  + theme_bw()    
poly1


## regress rii on gradient
m1 <- lm(Phacelia.biomass~arid.gradient, data=rii.mean)
summary(m1) ## not significant
## 
## Call:
## lm(formula = Phacelia.biomass ~ arid.gradient, data = rii.mean)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.022533 -0.013199 -0.005837  0.008594  0.058681 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.038659   0.020960   1.844   0.0899 .
## arid.gradient 0.008913   0.006828   1.305   0.2162  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02266 on 12 degrees of freedom
## Multiple R-squared:  0.1244, Adjusted R-squared:  0.05139 
## F-statistic: 1.704 on 1 and 12 DF,  p-value: 0.2162
m2 <- lm(Plantago.biomass~arid.gradient, data=rii.mean)
summary(m2) ## not significant
## 
## Call:
## lm(formula = Plantago.biomass ~ arid.gradient, data = rii.mean)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030560 -0.002438  0.002575  0.006742  0.019216 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.006458   0.014218  -0.454    0.658
## arid.gradient -0.001721   0.004632  -0.372    0.717
## 
## Residual standard error: 0.01537 on 12 degrees of freedom
## Multiple R-squared:  0.01137,    Adjusted R-squared:  -0.07101 
## F-statistic: 0.1381 on 1 and 12 DF,  p-value: 0.7167
m3 <- lm(Salvia.biomass~arid.gradient, data=rii.mean)
summary(m3) ## not significant
## 
## Call:
## lm(formula = Salvia.biomass ~ arid.gradient, data = rii.mean)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.051586 -0.015918 -0.001203  0.007707  0.066895 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.015346   0.029966  -0.512    0.618
## arid.gradient -0.004097   0.009762  -0.420    0.682
## 
## Residual standard error: 0.0324 on 12 degrees of freedom
## Multiple R-squared:  0.01447,    Adjusted R-squared:  -0.06766 
## F-statistic: 0.1761 on 1 and 12 DF,  p-value: 0.6821
# 
# m1 <- lm(Phacelia.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
# m2 <- lm(Plantago.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
# m3 <- lm(Salvia.biomass~arid.gradient, data=subset(rii.mean,Year == 2016))
# 
# 
# m1 <- lm(Phacelia.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
# m2 <- lm(Plantago.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
# m3 <- lm(Salvia.biomass~arid.gradient, data=subset(rii.mean,Year == 2017))
# 
# 
# plot(rii.mean[,"min.temp"], rii.mean[,"Phacelia.biomass"], ylim=c(-0.4,0.4), pch=19)
# points(rii.mean[,"min.temp"], rii.mean[,"Plantago.biomass"], pch=21, bg="blue")
# points(rii.mean[,"min.temp"], rii.mean[,"Salvia.biomass"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
# 
# rii.test <- subset(rii.mean, Gradient.x>3 & Year == 2017 | Year == 2016)
# 
# 
# plot(rii.test[,"arid.gradient"], rii.test[,"Phacelia.biomass"], ylim=c(-0.4,0.4), pch=19)
# points(rii.test[,"arid.gradient"], rii.test[,"Plantago.biomass"], pch=21, bg="blue")
# points(rii.test[,"arid.gradient"], rii.test[,"Salvia.biomass"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
# 
# 
# plot(rii.mean[,"arid.gradient"], rii.mean[,"Phacelia"], ylim=c(-0.4,0.4), pch=19)
# points(rii.mean[,"arid.gradient"], rii.mean[,"Plantago"], pch=21, bg="blue")
# points(rii.mean[,"arid.gradient"], rii.mean[,"Salvia"], pch=21, bg="red")
# abline(h=0, lty=2, lwd=2)
# 
# 
# m1 <- lm(Phacelia~ poly(arid.gradient,2), data=rii.mean)
# summary(m1)
# m2 <- lm(Plantago.biomass~arid.gradient, data=rii.mean)
# m3 <- lm(Salvia.biomass~arid.gradient, data=rii.mean)

Differences from Nutrient content

## join aridity with nutrients
nutrients.climate <- merge(nutrients,subset(site.climate, Year==2016), by=c("Site"))


## Nitrogen difference between shrub and sites
m.nit <- aov(log(N) ~ Site * microsite, data=nutrients)
summary(m.nit)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6  63.04   10.51  13.241 2.86e-09 ***
## microsite       1  40.16   40.16  50.614 2.24e-09 ***
## Site:microsite  6   4.96    0.83   1.042    0.408    
## Residuals      56  44.43    0.79                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.nit, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(N) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                diff      lwr      upr p adj
## shrub-open 1.514873 1.088317 1.941429     0
## Potassium difference between shrub and sites
m.pot <- aov(log(K) ~ Site * microsite, data=nutrients)
summary(m.pot)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 15.897   2.650  22.013 4.09e-13 ***
## microsite       1  4.445   4.445  36.934 1.14e-07 ***
## Site:microsite  6  1.605   0.268   2.223   0.0541 .  
## Residuals      56  6.740   0.120                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pot, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(K) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr       upr p adj
## shrub-open 0.5040081 0.3378755 0.6701406 1e-07
## Phosphorus difference between shrub and sites
m.pho <- aov(log(P) ~ Site * microsite, data=nutrients)
summary(m.pho)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## Site            6 23.163   3.861   14.33 8.10e-10 ***
## microsite       1 10.546  10.546   39.15 5.79e-08 ***
## Site:microsite  6  3.394   0.566    2.10   0.0676 .  
## Residuals      56 15.085   0.269                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m.pho, "microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(P) ~ Site * microsite, data = nutrients)
## 
## $microsite
##                 diff       lwr      upr p adj
## shrub-open 0.7762734 0.5277317 1.024815 1e-07
nutrient.mean <- nutrients.climate %>% group_by(Site, arid.gradient) %>% summarize(nitrogen=mean(N), phosphorus=mean(P), potassium=mean(K))
nutrient.mean <- data.frame(nutrient.mean)
 

ggplot(nutrient.mean) + geom_jitter(aes(x=arid.gradient, y=nitrogen), color="#E69F00", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=nitrogen), color="#E69F00") + geom_jitter(aes(x=arid.gradient, y=potassium/20), color="#56B4E9", size=3) +  stat_smooth(method="lm", formula= y~poly(x,2),aes(x=arid.gradient, y=potassium/20), color="#56B4E9") + geom_jitter(aes(x=arid.gradient, y=phosphorus), color="#009E73", size=3) +ylab("soil nutrient content")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"),panel.background = element_blank(),axis.text=element_text(size=14),axis.title=element_text(size=16)) + xlab("aridity gradient")+
annotate("segment", x = -4.5, xend = -4.35, y = 25, yend = 25,  colour = "#E69F00", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 25,  label = "Nitrogen", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 23, yend = 23,  colour = "#56B4E9", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 23,  label = "Potassium", size=5, hjust=0)+
annotate("segment", x = -4.5, xend = -4.35, y = 21, yend = 21,  colour = "#009E73", size=2) + ## add custom legend
  annotate("text", x = -4.3, y = 21,  label = "Phosphorus", size=5, hjust=0)

m1 <- lm(nitrogen~poly(arid.gradient,2), data=nutrient.mean)
summary(m1)
## 
## Call:
## lm(formula = nitrogen ~ poly(arid.gradient, 2), data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  1.1092  1.0570 -0.4565 -0.8511 -0.8789 -0.3757  0.3960 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               5.4326     0.3942  13.782 0.000161 ***
## poly(arid.gradient, 2)1 -12.6364     1.0429 -12.116 0.000266 ***
## poly(arid.gradient, 2)2   7.0832     1.0429   6.792 0.002454 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.043 on 4 degrees of freedom
## Multiple R-squared:  0.9797, Adjusted R-squared:  0.9695 
## F-statistic: 96.47 on 2 and 4 DF,  p-value: 0.0004126
m2 <- lm(phosphorus~arid.gradient, data=nutrient.mean)
summary(m2)
## 
## Call:
## lm(formula = phosphorus ~ arid.gradient, data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  4.3035  1.1485 -2.0485  1.1930 -0.2791 -5.5020  1.1845 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     14.191      4.527   3.135   0.0258 *
## arid.gradient    1.780      1.298   1.371   0.2287  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.383 on 5 degrees of freedom
## Multiple R-squared:  0.2732, Adjusted R-squared:  0.1278 
## F-statistic: 1.879 on 1 and 5 DF,  p-value: 0.2287
m3 <- lm(potassium~poly(arid.gradient,2), data=nutrient.mean)
summary(m3)
## 
## Call:
## lm(formula = potassium ~ poly(arid.gradient, 2), data = nutrient.mean)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  -3.551 -65.927  32.611  57.492 -25.220  52.524 -47.931 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)   
## (Intercept)               186.56      22.69   8.221  0.00119 **
## poly(arid.gradient, 2)1   189.58      60.04   3.157  0.03427 * 
## poly(arid.gradient, 2)2   190.31      60.04   3.170  0.03387 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 60.04 on 4 degrees of freedom
## Multiple R-squared:  0.8334, Adjusted R-squared:  0.7502 
## F-statistic: 10.01 on 2 and 4 DF,  p-value: 0.02774

annual plant community repsonse

spp.data  <- read.csv("Data/ERG.communitydata.csv")
spp.data[is.na(spp.data)] <- 0

mean.spp <- spp.data %>% group_by(Year,Site, Microsite) %>%  summarize(abd=mean(Abundance), richness=mean(Richness), biomass=mean(Biomass))
mean.spp <- data.frame(mean.spp)

## collect aridity gradient data
site.climate <- rbind(s1.mean,s2.mean)
site.climate[,"Year"] <- c(rep("2016",7),rep("2017",7))

## combine climate data with community data
mean.spp <- merge(mean.spp, site.climate, by.y=c("Year","Site"))

# ## community responses season1
# mean.spp2016 <- subset(mean.spp, Year==2016)
# 
# ## richness
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","richness"], ylab="average site richness", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","richness"], pch=19) 
# 
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2016)
# summary(m1)
# 
# ## abundance
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","abd"], pch=19) 
# 
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2016)
# summary(m2)
# 
# ## biomass
# plot(mean.spp2016[mean.spp2016$Microsite=="open","arid.gradient"], mean.spp2016[mean.spp2016$Microsite=="open","biomass"], ylab="average site biomass", xlab="aridity gradient")
# points(mean.spp2016[mean.spp2016$Microsite=="shrub","arid.gradient"],mean.spp2016[mean.spp2016$Microsite=="shrub","biomass"], pch=19) 
# 
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2016)
# summary(m3)
# 
# ##season2
# 
# ## community responses season2
# 
# mean.spp2017 <- subset(mean.spp, Year==2017)
# 
# 
# ## richness
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","richness"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","richness"], pch=19) 
# 
# m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp2017)
# summary(m1)
# 
# ## abundance
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","abd"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","abd"], pch=19) 
# 
# m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp2017)
# summary(m2)
# 
# ## biomass
# plot(mean.spp2017[mean.spp2017$Microsite=="open","arid.gradient"], mean.spp2017[mean.spp2017$Microsite=="open","biomass"])
# points(mean.spp2017[mean.spp2017$Microsite=="shrub","arid.gradient"],mean.spp2017[mean.spp2017$Microsite=="shrub","biomass"], pch=19) 
# 
# m3 <- lm(biomass~arid.gradient*Microsite, data=mean.spp2017)
# summary(m3)


## both seasons

## richness
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","richness"],  ylab="average site richness", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","richness"], pch=19) 

m1 <- lm(richness~arid.gradient*Microsite, data=mean.spp)
summary(m1)
## 
## Call:
## lm(formula = richness ~ arid.gradient * Microsite, data = mean.spp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56169 -0.66225 -0.07639  0.66154  1.76433 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    4.9206     0.8381   5.871 4.68e-06 ***
## arid.gradient                  0.8699     0.2730   3.186  0.00397 ** 
## Micrositeshrub                -1.8968     1.1852  -1.600  0.12260    
## arid.gradient:Micrositeshrub  -0.5531     0.3861  -1.432  0.16490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9061 on 24 degrees of freedom
## Multiple R-squared:  0.3356, Adjusted R-squared:  0.2526 
## F-statistic: 4.042 on 3 and 24 DF,  p-value: 0.01849
## abundance
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], mean.spp[mean.spp$Microsite=="open","abd"], ylab="average site abundance", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],mean.spp[mean.spp$Microsite=="shrub","abd"], pch=19) 

m2 <- lm(abd~arid.gradient*Microsite, data=mean.spp)
summary(m2)
## 
## Call:
## lm(formula = abd ~ arid.gradient * Microsite, data = mean.spp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.997 -26.611  -5.081  25.895 108.541 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    162.65      38.28   4.249  0.00028 ***
## arid.gradient                   41.28      12.47   3.310  0.00294 ** 
## Micrositeshrub                  39.22      54.14   0.724  0.47582    
## arid.gradient:Micrositeshrub    10.04      17.64   0.569  0.57457    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.39 on 24 degrees of freedom
## Multiple R-squared:  0.5409, Adjusted R-squared:  0.4835 
## F-statistic: 9.426 on 3 and 24 DF,  p-value: 0.000268
## biomass
plot(mean.spp[mean.spp$Microsite=="open","arid.gradient"], ihs(mean.spp[mean.spp$Microsite=="open","biomass"]), ylab="average site biomass", xlab="aridity gradient")
points(mean.spp[mean.spp$Microsite=="shrub","arid.gradient"],ihs(mean.spp[mean.spp$Microsite=="shrub","biomass"]), pch=19)

m3 <- lm(ihs(biomass)~arid.gradient*Microsite, data=mean.spp)
summary(m3)
## 
## Call:
## lm(formula = ihs(biomass) ~ arid.gradient * Microsite, data = mean.spp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6596 -0.5085 -0.0684  0.3314  1.2077 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    2.6251     0.5577   4.707 8.75e-05 ***
## arid.gradient                  0.5873     0.1817   3.233  0.00355 ** 
## Micrositeshrub                 1.6603     0.7888   2.105  0.04596 *  
## arid.gradient:Micrositeshrub   0.3698     0.2569   1.439  0.16303    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.603 on 24 degrees of freedom
## Multiple R-squared:  0.6498, Adjusted R-squared:  0.606 
## F-statistic: 14.84 on 3 and 24 DF,  p-value: 1.128e-05

cleaned up community response to gradient

RDA community and environmental data

community <- read.csv("Data/ERG.communitydata.csv")
community[is.na(community)] <- 0

##2016
## add enviro data to dataframe
comm2016 <- subset(community, Year==2016)
site.vars2016[,"Site"] <- as.factor(row.names(site.vars2016))
comm2016 <- merge(comm2016, site.vars2016, by="Site")

## transform spp data
community.trans <- decostand(comm2016[,13:53], "hellinger")

rda1 <- rda(community.trans, comm2016[,54:57])
(R2adj <- RsquareAdj(rda1)$adj.r.squared) ## 50.4% of variation explained
## [1] 0.3871872
## build byplot manually
par(mar=c(4.5,4.5,0.5,0.5))
plot(rda1, type="n", xlim=c(-1,1), ylim=c(-1.5,1.5))

## calculate priority
spp.priority <- colSums(comm2016[,13:53])

## plot RDA1
colvec <- rep(c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2"),each=60)
points(rda1, display = "sites", col = "black", pch = c(21), bg = colvec)
orditorp(rda1, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=0.5)
text(rda1,  display = "bp", col = "blue", cex = 0.8) 
legend("bottomright", pch=c(21), legend=unique(comm2016$Site), pt.bg=unique(colvec), cex=1)

## ordikplot to adjust species locations

## 2017
## add enviro data to dataframe
comm2017 <- subset(community, Year==2017)
site.vars[,"Site"] <- as.factor(row.names(site.vars))
comm2017 <- merge(comm2017, site.vars, by="Site")
# 
# ## transform spp data
# community.trans <- decostand(comm2017[,13:53], "hellinger")
# 
# rda2 <- rda(community.trans, comm2017[,54:58])
# (R2adj <- RsquareAdj(rda2)$adj.r.squared) ## 50.4% of variation explained
# 
# 
# ## build byplot manually
# par(mar=c(4.5,4.5,0.5,0.5))
# plot(rda2, type="n", xlim=c(-1,1), ylim=c(-1.5,1.5))
# 
# ## calculate priority
# spp.priority <- colSums(comm2017[,13:53])
# 
# ## plot RDA2
# colvec <- rep(c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2"),each=60)
# points(rda2, display = "sites", col = "black", pch = c(21), bg = colvec)
# orditorp(rda2, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=0.5)
# text(rda2,  display = "bp", col = "blue", cex = 0.8) 
# legend("bottomright", pch=c(21), legend=unique(comm2016$Site), pt.bg=unique(colvec), cex=1)
# ## ordikplot to adjust species locations

Phytometer analysis

census <- read.csv("Data/ERG.phytometer.census.csv")
census[is.na(census)] <- 0
census[,"phyto.abd"] <- rowSums(census[,c("Phacelia","Plantago","Salvia")])

## calculate the number of seeds per gram
seed.mass <- read.csv("Data/seed.mass.csv")
seed.mass.avg <- seed.mass %>% group_by(species) %>%  summarize(avg.seed.gram=mean(seed.number),se.seed.gram=se(seed.number))
seed.mass.avg <- data.frame(seed.mass.avg)
seed.mass.avg
##      species avg.seed.gram se.seed.gram
## 1   Amsinkia         235.2     7.611833
## 2 Caulanthus        3019.0     8.916277
## 3  Lipidium          752.4     4.261455
## 4 Monolopia          737.8     7.831986
## 5   Phacelia         772.0     5.300943
## 6   Plantago         558.2     6.865858
## 7     Salvia         980.0     9.612492
## get proportion of seed germinated per 0.3 grams of seed
census[,"Phacelia.prop"] <- census[,"Phacelia"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Phacelia")]*0.3)
census[,"Plantago.prop"] <- census[,"Plantago"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Plantago")]*0.3)
census[,"Salvia.prop"] <- census[,"Salvia"]/(seed.mass.avg$avg.seed.gram[which(seed.mass.avg$species=="Salvia")]*0.3)


## beginning
census.int <- subset(census, Census=="emergence")

## all species
m1 <- glm.nb(Phacelia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m1, test="Chisq") ## Site * Micro and Year significant
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.8108), link: log
## 
## Response: Phacelia
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                      839    1079.71              
## Microsite                1   10.304       838    1069.41  0.001328 ** 
## Nutrient                 1    1.234       837    1068.17  0.266538    
## Site                     6  159.651       831     908.52 < 2.2e-16 ***
## Year                     1   31.424       830     877.10 2.073e-08 ***
## Microsite:Nutrient       1    0.046       829     877.05  0.829874    
## Microsite:Site           6   21.032       823     856.02  0.001811 ** 
## Nutrient:Site            6    5.542       817     850.48  0.476389    
## Microsite:Nutrient:Site  6    8.030       811     842.45  0.235879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Phacelia
m1 <- glm.nb(Phacelia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m1, test="Chisq") ## Site * Micro and Year significant
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.8108), link: log
## 
## Response: Phacelia
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                      839    1079.71              
## Microsite                1   10.304       838    1069.41  0.001328 ** 
## Nutrient                 1    1.234       837    1068.17  0.266538    
## Site                     6  159.651       831     908.52 < 2.2e-16 ***
## Year                     1   31.424       830     877.10 2.073e-08 ***
## Microsite:Nutrient       1    0.046       829     877.05  0.829874    
## Microsite:Site           6   21.032       823     856.02  0.001811 ** 
## Nutrient:Site            6    5.542       817     850.48  0.476389    
## Microsite:Nutrient:Site  6    8.030       811     842.45  0.235879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Plantago
m2 <- glm.nb(Plantago~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m2, test="Chisq") ## Site * Micro and Year significant
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.5335), link: log
## 
## Response: Plantago
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                      839     984.30              
## Microsite                1    6.256       838     978.04   0.01238 *  
## Nutrient                 1    1.060       837     976.98   0.30332    
## Site                     6  162.687       831     814.30 < 2.2e-16 ***
## Year                     1    4.138       830     810.16   0.04194 *  
## Microsite:Nutrient       1    2.388       829     807.77   0.12224    
## Microsite:Site           6   43.000       823     764.77 1.167e-07 ***
## Nutrient:Site            6    6.391       817     758.38   0.38088    
## Microsite:Nutrient:Site  6    4.615       811     753.77   0.59403    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Salvia
m3 <- glm.nb(Salvia~ Microsite * Nutrient * Site + Year, data=census.int)
anova(m3, test="Chisq") ## Site * Micro and Year significant + (micro x nutrient)
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.8273), link: log
## 
## Response: Salvia
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                      839    1210.25              
## Microsite                1    5.924       838    1204.33 0.0149386 *  
## Nutrient                 1    1.513       837    1202.82 0.2186678    
## Site                     6  211.625       831     991.19 < 2.2e-16 ***
## Year                     1   17.463       830     973.73 2.929e-05 ***
## Microsite:Nutrient       1    7.824       829     965.90 0.0051548 ** 
## Microsite:Site           6   26.421       823     939.48 0.0001858 ***
## Nutrient:Site            6    4.587       817     934.90 0.5977624    
## Microsite:Nutrient:Site  6    1.566       811     933.33 0.9549811    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
census.int.plot <- gather(census.int, species, abundance, Phacelia:Salvia)

##calculate confidence interval
census.plot<- census.int.plot %>%  group_by(Microsite, species, Site, Year) %>%  summarize(avg=mean(abundance),ci=se(abundance)*1.96)

ggplot(census.plot, aes(x=Microsite, y=avg, fill=species))+
  geom_bar(position=position_dodge(), stat="identity")+
   geom_errorbar(aes(ymin=avg-ci, ymax=avg+ci),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ scale_fill_brewer() +ylab("Abundance")+ theme_Publication() + facet_grid(Year~Site) 

## end season
census.end <- subset(census, Census=="end")
census.end[,"Year"] <- as.factor(census.end$Year)

## run full model for phyto biomass
m1 <- aov(log(phyto.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, phyto.biomass>0))
summary(m1)
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6  330.8   55.13  32.164  < 2e-16 ***
## Microsite                 1   10.9   10.88   6.347 0.012127 *  
## Nutrient                  1   30.3   30.25  17.651 3.24e-05 ***
## Year                      1  161.2  161.21  94.051  < 2e-16 ***
## Site:Microsite            6   23.9    3.99   2.327 0.031940 *  
## Site:Nutrient             6   41.9    6.98   4.073 0.000551 ***
## Microsite:Nutrient        1    1.7    1.73   1.011 0.315196    
## Site:Microsite:Nutrient   6    7.9    1.32   0.768 0.595068    
## Residuals               424  726.7    1.71                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m1, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(phyto.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, phyto.biomass > 0))
## 
## $Microsite
##                 diff        lwr       upr     p adj
## shrub-open 0.3096905 0.06750868 0.5518722 0.0123243
## all species end of season
m2 <- glm.nb(phyto.abd~ Microsite * Nutrient * Site + Year, data=census.end)
anova(m2, test="Chisq") ## Site and microsite*nutrient significant
## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.455), link: log
## 
## Response: phyto.abd
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                                      839    1009.51             
## Microsite                1    1.059       838    1008.45  0.30341    
## Nutrient                 1    2.527       837    1005.92  0.11194    
## Site                     6  135.420       831     870.50  < 2e-16 ***
## Year                     1    3.412       830     867.09  0.06474 .  
## Microsite:Nutrient       1    4.586       829     862.50  0.03224 *  
## Microsite:Site           6    5.340       823     857.16  0.50104    
## Nutrient:Site            6    5.169       817     852.00  0.52238    
## Microsite:Nutrient:Site  6    5.683       811     846.31  0.45956    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## run full model for phacelia biomass
m3 <- aov(log(Phacelia.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Phacelia.biomass>0))
summary(m3) ## site and microsite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6   89.7   14.95   8.274 4.95e-08 ***
## Microsite                 1   21.2   21.20  11.738  0.00074 ***
## Nutrient                  1   10.2   10.19   5.642  0.01845 *  
## Year                      1   95.5   95.49  52.857 7.42e-12 ***
## Site:Microsite            6   23.6    3.93   2.173  0.04692 *  
## Site:Nutrient             6    8.0    1.34   0.742  0.61612    
## Microsite:Nutrient        1    0.0    0.00   0.001  0.97552    
## Site:Microsite:Nutrient   5    4.1    0.82   0.452  0.81122    
## Residuals               205  370.3    1.81                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m3, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Phacelia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Phacelia.biomass > 0))
## 
## $Microsite
##                 diff       lwr       upr     p adj
## shrub-open 0.6312547 0.2667053 0.9958041 0.0007714
## run full model for plantago biomass
m4 <- aov(log(Plantago.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Plantago.biomass>0))
summary(m4) ## site and year + sitexmicrosite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6 116.35   19.39  19.843 2.32e-16 ***
## Microsite                 1   1.96    1.96   2.001    0.160    
## Nutrient                  1  32.03   32.03  32.778 6.80e-08 ***
## Year                      1 138.38  138.38 141.597  < 2e-16 ***
## Site:Microsite            6   4.87    0.81   0.830    0.549    
## Site:Nutrient             4   1.99    0.50   0.508    0.730    
## Microsite:Nutrient        1   3.17    3.17   3.244    0.074 .  
## Site:Microsite:Nutrient   3   0.65    0.22   0.220    0.882    
## Residuals               130 127.05    0.98                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m4, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Plantago.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Plantago.biomass > 0))
## 
## $Microsite
##                  diff        lwr       upr     p adj
## shrub-open -0.2179629 -0.5348829 0.0989571 0.1759823
## run full model for plantago biomass
m5 <- aov(log(Salvia.biomass)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Salvia.biomass>0))
summary(m5) ## site and year + sitexmicrosite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6  170.1   28.35  19.982  < 2e-16 ***
## Microsite                 1    0.0    0.01   0.008 0.929221    
## Nutrient                  1   21.4   21.43  15.103 0.000124 ***
## Year                      1   44.1   44.05  31.045 5.39e-08 ***
## Site:Microsite            6    6.8    1.14   0.802 0.568988    
## Site:Nutrient             6   46.2    7.70   5.424 2.35e-05 ***
## Microsite:Nutrient        1    0.5    0.55   0.387 0.534459    
## Site:Microsite:Nutrient   6   12.9    2.16   1.520 0.171069    
## Residuals               317  449.8    1.42                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m5, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Salvia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Salvia.biomass > 0))
## 
## $Microsite
##                  diff        lwr       upr     p adj
## shrub-open 0.01133515 -0.2406604 0.2633307 0.9295351
## run full model for phacelia abundance
m6 <- glm.nb(Phacelia.biomass~ Site * Microsite * Nutrient + Year, data=subset(census.end, Phacelia>0))
summary(m3) ## site and microsite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6   89.7   14.95   8.274 4.95e-08 ***
## Microsite                 1   21.2   21.20  11.738  0.00074 ***
## Nutrient                  1   10.2   10.19   5.642  0.01845 *  
## Year                      1   95.5   95.49  52.857 7.42e-12 ***
## Site:Microsite            6   23.6    3.93   2.173  0.04692 *  
## Site:Nutrient             6    8.0    1.34   0.742  0.61612    
## Microsite:Nutrient        1    0.0    0.00   0.001  0.97552    
## Site:Microsite:Nutrient   5    4.1    0.82   0.452  0.81122    
## Residuals               205  370.3    1.81                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m3, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Phacelia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Phacelia.biomass > 0))
## 
## $Microsite
##                 diff       lwr       upr     p adj
## shrub-open 0.6312547 0.2667053 0.9958041 0.0007714
## run full model for plantago biomass
m7 <- aov(log(Plantago)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Plantago>0))
summary(m4) ## site and year + sitexmicrosite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6 116.35   19.39  19.843 2.32e-16 ***
## Microsite                 1   1.96    1.96   2.001    0.160    
## Nutrient                  1  32.03   32.03  32.778 6.80e-08 ***
## Year                      1 138.38  138.38 141.597  < 2e-16 ***
## Site:Microsite            6   4.87    0.81   0.830    0.549    
## Site:Nutrient             4   1.99    0.50   0.508    0.730    
## Microsite:Nutrient        1   3.17    3.17   3.244    0.074 .  
## Site:Microsite:Nutrient   3   0.65    0.22   0.220    0.882    
## Residuals               130 127.05    0.98                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m4, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Plantago.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Plantago.biomass > 0))
## 
## $Microsite
##                  diff        lwr       upr     p adj
## shrub-open -0.2179629 -0.5348829 0.0989571 0.1759823
## run full model for plantago biomass
m8 <- aov(log(Salvia)~ Site * Microsite * Nutrient + Year, data=subset(census.end, Salvia>0))
summary(m5) ## site and year + sitexmicrosite
##                          Df Sum Sq Mean Sq F value   Pr(>F)    
## Site                      6  170.1   28.35  19.982  < 2e-16 ***
## Microsite                 1    0.0    0.01   0.008 0.929221    
## Nutrient                  1   21.4   21.43  15.103 0.000124 ***
## Year                      1   44.1   44.05  31.045 5.39e-08 ***
## Site:Microsite            6    6.8    1.14   0.802 0.568988    
## Site:Nutrient             6   46.2    7.70   5.424 2.35e-05 ***
## Microsite:Nutrient        1    0.5    0.55   0.387 0.534459    
## Site:Microsite:Nutrient   6   12.9    2.16   1.520 0.171069    
## Residuals               317  449.8    1.42                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(m5, "Microsite")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = log(Salvia.biomass) ~ Site * Microsite * Nutrient + Year, data = subset(census.end, Salvia.biomass > 0))
## 
## $Microsite
##                  diff        lwr       upr     p adj
## shrub-open 0.01133515 -0.2406604 0.2633307 0.9295351
census.end.plot <- gather(census.end, species, biomass, Phacelia.biomass:Salvia.biomass)

##calculate confidence interval
census.plot<- census.end.plot %>%  group_by(Microsite, species, Site, Year) %>%  summarize(avg=mean(biomass),ci=se(biomass)*1.96)

ggplot(census.plot, aes(x=Microsite, y=avg, fill=species))+
  geom_bar(position=position_dodge(), stat="identity")+
   geom_errorbar(aes(ymin=avg-ci, ymax=avg+ci),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+ scale_fill_brewer() +ylab("Abundance")+ theme_Publication() + facet_grid(Year~Site) 

### RDA
# 
# ##collect environmental variables for site
# nutrients <- read.csv("Data/ERG.soilnutrients.csv")
# nutrients.mean <- nutrients %>% group_by(gradient,site, microsite) %>% summarise_each(funs(mean))
# nutrients.vars <- data.frame(nutrients.mean)
# 
# ## minimum explanation of variables
# ambient2016[is.na(ambient2016)] <- 0
# community <- ambient2016 %>% group_by(Gradient,Site, Microsite) %>% summarise_each(funs(sum))
# community <- data.frame(community)
# 
# end.census <- subset(census2016, census=="end")
# end.census  <- end.census  %>% group_by(Gradient,Site, Microsite) %>% summarise_each(funs(mean))
# end.census<- data.frame(end.census)
# 
# ## daily means
# day.mean <- aggregate(HOBO.data, by=list(day=HOBO.data$Day, micro=HOBO.data$Microsite, site=HOBO.data$Site), mean)
# 
# ## summarize daily means across sites
# means <- day.mean %>% group_by(gradient, micro) %>% summarize(temp=mean(Temp),rh=mean(RH),temp.se=se(Temp),rh.se=se(RH))
# means <- data.frame(means)
# 
# 
# envs <- data.frame(swc=end.census[,"swc"],nutrients.vars[,c("N","P","K")], means[,c("temp","rh")])
# 
# ## hellinger transformation
# community.trans <- decostand(community[,12:42], "hellinger")
# 
# ## rda with environmental variables
# rda1 <- rda(community.trans, envs)
# anova(rda1)
# (R2adj <- RsquareAdj(rda1)$adj.r.squared) ## 25.7% of variation explained
# 
# 
# ## build byplot manually
# par(mar=c(4.5,4.5,0.5,0.5))
# plot(rda1, type="n", xlim=c(-1,1))
# 
# with(community, levels(Site))
# 
# spp.priority <- colSums(community[,12:42])
# 
# colvec <- c("blue","dodgerblue3","cyan","yellow2","orange","tomato","red2")
# with(community, points(rda1, display = "sites", col = "black", pch = c(21,22), bg = colvec[Site]))
# orditorp(rda1, display = "species", cex = 0.7, col = "darkred", priority=spp.priority, air=1,)
# text(rda1,  display = "bp", col = "blue", cex = 0.8) 
# legend("bottomright", pch=c(21), legend=community$Site[community$Microsite=="shrub"], pt.bg=colvec, cex=1)
# 
# ## check variable explanation
# vif.cca(rda1)
# 
# v.clim <- cbind(envs[,c("temp","rh")])
# v.nut <- cbind(envs[,c("N","P","K")])
# v.swc <- envs[,"swc"]
# 
# x <- varpart(community.trans,v.clim,v.nut,v.swc)
# 
# showvarparts(3)

Difference in effect size

library(bootES)
## Loading required package: boot
## 
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
## 
##     melanoma
effect.cal <- function(x){ bootES(x, R=999, data.col="Biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}

effect.site <- by(subset(comm2016, Site != "Barstow"), comm2016$Site[comm2016$Site != "Barstow"], FUN=effect.cal)

site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa)) 
site.summary[5,] <- 0
site.summary <- sapply(site.summary, as.numeric)

par(mar=c(4.5,4.5,.5,.5))
plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
arrows(c(1:7),as.numeric(site.summary[,"ci.high"]),c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
abline(h=0, lty=2, lwd=2)
axis(1, c(1:7), labels=site.names, cex.axis=1)

## phytometer

effect.cal <- function(x){ bootES(x, R=999, data.col="phyto.biomass", group.col="Microsite", contrast=c(shrub=1,open=-1), effect.type=c("cohens.d"))}

census2016 <- subset(census, Census=="end" & Year ==2016)

effect.site <- by(subset(census2016, Site != "Barstow"), census2016$Site[census2016$Site != "Barstow"], FUN=effect.cal)

site.summary <- rbind(summary(effect.site$Panoche),summary(effect.site$Cuyama),summary(effect.site$TejonRanch),data.frame(stat=0,ci.low=0,ci.high=0,bias=0,std.error=0),summary(effect.site$MojavePreserve),summary(effect.site$SheepholeValley),summary(effect.site$Tecopa)) 
site.summary <- sapply(site.summary, as.numeric)

par(mar=c(4.5,4.5,.5,.5))
plot(c(1:7),site.summary[,"stat"], xlim=c(0.5,7.5), ylim=c(-1,3), pch=19, cex=1.4, ylab="Hedge's G", xlab="", xaxt="n")
arrows(c(1:7),site.summary[,"ci.high"],c(1:7),site.summary[,"ci.low"], angle=90, code=3, length=0, lwd=2)
abline(h=0, lty=2, lwd=2)
axis(1, c(1:7), labels=site.names, cex.axis=1)

Rii between years

rii.dat <- rii(community, j=1:8, var=c("Abundance","Richness","Biomass"))

rii.mean <- rii.dat %>% group_by(Year, Gradient, Site) %>% summarize(avg=mean(Biomass),se=se(Biomass))
rii.mean <- data.frame(rii.mean)

plot(rii.mean[rii.mean$Year==2016,"Gradient"]-.1,rii.mean[rii.mean$Year==2016,"avg"], ylim=c(-1,1), pch=19, cex=1.5, xlim=c(0.5,7.5), ylab="Rii Biomass", xlab="Site", xaxt="n", cex.lab=1.5)
axis(1, 1:7, lab=unique(rii.mean$Site))
error.bar(rii.mean[rii.mean$Year==2016,"Gradient"]-.1,rii.mean[rii.mean$Year==2016,"avg"],rii.mean[rii.mean$Year==2016,"se"])
error.bar(rii.mean[rii.mean$Year==2017,"Gradient"]+.1,rii.mean[rii.mean$Year==2017,"avg"],rii.mean[rii.mean$Year==2017,"se"])
points(rii.mean[rii.mean$Year==2017,"Gradient"]+.1,rii.mean[rii.mean$Year==2017,"avg"], pch=21, bg="Grey50", cex=1.5)
abline(h=0, lwd=2, lty=2)
legend(6.3,-0.55, c("2016","2017"), pch=22, pt.bg=c("Black","Grey50"), cex=1.6)

## add year column to precipitation data
precip[,"Year"] <- ifelse(precip$season=="season.1","2016","2017")
rii.precip <- merge(precip,rii.mean, by=c("Year","Gradient"))